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ABSTRACT 


The main objective of this project is to introduce the theory of queueing systems 
and to demonstrate its applicability to real life problems. In the first chapter, we discuss the 
Markovian property and measures of effectiveness for queueing systems with exponential 
interarrival and service times. This information is used to optimize the performance of 
the given system. The second chapter of this project is concerned with queueing systems 
with exponential interarrival times, Erlang service times, and a single server. The system 
performance measures for this queueing system are derived. An appropriate example model 
is constructed and investigated. Chapter three discusses different goodness-of-fit tests that 
can be used to determine whether the exponential distribution is appropriate for a given set 
of data. Advantages and disadvantages of using different goodness-of-fit tests are.discussed 
and relevant examples are given. Simulation is a method of analyzing queueing systems 
which is an alternative to the analytical approach. In the fourth chapter, a single server 
queueing system with exponential interarrival times and Erlang service times is simulated 


using Visual Basic for Applications. 
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Chapter 1 


Introduction 


1.1 Overview 


The study of the characteristics of waiting lines, or queues, has many important 
applications. One of the first problems studied in the field of queueing theory was telephone 
traffic congestion by A. K. Erlang in 1909. Erlang’s research sparked interest among 
many other mathematicians who extended his work. Up to the 1950’s, telephony was 
the principal application of this theory. Today, queueing theory is a useful tool in other 
important applications such as air traffic control, machine repair, scheduling problems, and 


time-sharing. 


Terminology and Kendall-Lee Notation 


To describe a queueing system, we must specify the arrival and service processes. 
The arrival process of most queueing systems is independent of the number of customers 
present and may be described by the probability distribution that governs the time between 
successive arrivals, called the interarrival time. We will assume that at most one arrival 
may occur at any given instant and that the arrival process does not depend on the number 
of customers present in the queueing system. 

The service process of many systems is also independent of the number of cus- 
tomers present, and we may specify the service time distribution which describes a cus- 
tomer’s service time. We will assume that the servers do not work any slower or faster 


depending on the queue length. In this project, we will consider systems with a single 


queue and parallel servers. Servers are in parallel if each customer may be fully served by 
any server and then leave the system. 
In this work, we will use the following standard abbreviations for the probability 


distributions to describe interarrival and service times: | 


M = Times are independent, identically distributed random | 
variables having an exponential distribution. . 2 wa 
E, = Times are independent, identically distributed random 


variables having an Erlang distribution with shape parameter k. 


There is a standard notation created by Kendall (1951) for describing many queueing 
systems. This notation is used for models in which all arriving customers wait in a single 
queue until one of s parallel servers becomes available. There are six characteristics of the 
notation written in the form 1/2/3/ 4/5/6. The first characteristic describes the nature of 
the arrival process. The second describes the nature of the service times. | 
The third characteristic is the number, s, of parallel servers. The fourth char- 
_ acteristic describes the queue discipline. In this work, the order in which customers are 
served is first come, first served (FCF'S) . 
The fifth characteristic i is the maximum number of customers allowed in ies sys- 
tem, including customers in line and customers being served. If there is no maximum, we 
denote this by oo. The sixth characteristic is the size of the population from which cus- 
tomers are drawn. We denote this by co unless the number of potential customers is of the 
“same order of magnitude as the number of Leas servers. is for 1/2/s/FCFS/oo / oO, 
the shorthand notation of 1 /2/s is used. 


System Performance Parameters © 


Certain characteristics of a queueing system are of particular interest to optimize 
its performance. The most important characteristics are as follows: | 
1; = steady-state probability that j customers are in the system 
L = expected number of customers in the system 
Lg = expected number of customers in queue 


Ls = expected number of customers in service 


W= expected time a customer spends in 1 the system 
Wy= expected time. a customer spends i in queue 
W, = expected time a customer ‘spends i in service 
A= average number of customers per unit time (arrival rate) 
‘pb = average number of service. completions per unit time (service, rate) 


p = /sp = traffic intensity. 


1.2. Birth-Death Processes , 


We assume that no. more than one arrival 1 may occur at any given instant of time. 
Define t; as the time of arrival of the ith customer. Let T; = dea _ ti be the ith interarrival. 
time. We will assume the T;’s are independent, continuous random variables described by 
the random variable A. . 

‘The most appropriate distribution for A is the enonential distribution since inter- 
arrival times are not usually very long. The density function of an ee distribution 
with parameter d is : 
a(t) = de®, ai 


_ where 2 is the abival rate, measured i in units of arrivals per. chon Also, 1/2 is defined 
as the mean interarrival time. 2 
We know from [Ros02] that the mean and variance of an exponcnuals random 


variable are given by 


BCA) 


aes 


varA 


In this chapter, we assume 2 that service. times are exponential. Let Ki be the service ‘rate 


- in units of customers per hour. ‘Then, 5 y, pis the mean service time, 


Another reason that the exponential. distribution is often used for interarrival 
- times is because of its: Markovian property. This property implies that the probability 
distribution of the time until the next: arrival does not depend on how long it has been 
since the last arrival. Throughout | this work, we assume that. all interarrival times are 


. exponential. The proof of the following lemma can n be found i in [Win94. 


Lemma 1.1 (Markovian Property). If A has an exponential distribution, then for all 


nonnegative values of t and h 
P(A>t+h|A>t)=P(A>h). 
The proof of the following theorem can be found in [GH85]. 


Theorem 1.2. Consider an arrival process {N(t), t > 0}, where N(t) is the number of 
arrivals up to t, N(0) =0, and satisfies the following assumptions: 

(i) P(an arrival occurs between time t and t+ At) = AAt + o(At), where d is a constant 
independent of N(t) and lim o(At)/t = 0; 

(it) P(more than one arrival between t and t + At) = o(At); 

(iii) the number of arrivals in non-overlapping time intervals are statistically independent. 
Then, interarrival times are exponential with parameter X if the number of arrivals that 


occur in an interval of length t follows a Poisson distribution with parameter At. 


Define the number of people present in a queueing system at time ¢ to be the state 
of the queueing system at time t. The probability of 7 people present in a queueing system 
at time ¢ given i people initially present is denoted by P;; (¢). For many queueing systems, 
P;; (¢) will approach a limit 7; for large t. This limit is independent of the number of people 
initially in the system. We call 2; the steady-state probability of state j, or alternatively 
the equilibrium probability of state 7. We will think of 7; as the probability that at any 
given time in the distant future, j7 customers will be present in the system. The 7; may 
also be interpreted as the fraction of time that 7 customers are present in the queueing 
system for some time in the distant future. 

A birtl-death process satisfies the following three laws: 

(i) A birth occurs between time ¢ and t+ At with probability 4;At + o(Az), where A; is 
called the birth rate in state 7. Hence a birth increases the state of the system from 7 to 
g+1. A birth is equivalent to an arrival to the system. 

(ii) A death occurs between time ¢t and t+ At with probability 4;At+ o(At), where p; i 
called the death rate in state 7. Thus a death decreases the system’s state from 7 to 7-1. 
A death will be thought of as a service completion. To ensure a non-negative state, we 
require ji9 = 0. 


(iii) Births and deaths are independent of each other. 


The birth-death model is not appropriate when either the interarrival or service 


times are not exponential. 


Steady-State Probabilities for Birth-Death Processes 


We will outline the steps of the derivation of the steady-state probabilities. See 
[Win94] for more details. 

We relate P,j (t + At) to Py (¢) for small At by 
Pay (t+ At) = Py (@) + At Ay-1Pig—1 @) + Hj41 Pigs (t) — Pig (t) hy — Pig (4) Ag) +0 (At). 
We divide both sides by A?é, and take the limit as At — oo. For all i and j > 1, we have 


Pig (t) = Aga Pi,g—a (t) + oy 41Pij41 (1) — Pig () oy — Pig As, 
and for 7 = 0, 
Pio (t) = w1P,1 (t) — AoPio (t) - 
We will use this infinite system of differential equations to obtain the steady-state proba- 
bilities defined by 
T= jim Pi; (t). 
To find the steady-state, we set P,;’(t) = 0 and acquire 
mit THM = ™HAYtHs) G2) 


T1H1 = ToAo. 


From the above equations, it follows that 


Too 

1 o> 

M1 
ee 70 (AoA1) 
Hip 

Let 
7 AoA + Aj—1 
Cj = 

Biba: :> Bj 


Then, 


Tj = Toc;. 1.1 
w) J 


In the next section, we apply the theory of birth-death processes to determine 
the steady-state probabilities for the M/M/s queueing system and use them to derive 


measures of effectiveness for the system. 


1.3 Optimization for the M/M/s System 


Analysis of queueing systems may be used to answer one of the most important 
questions for the system: How can we minimize costs for a particular queueing model? 
Also, what is the optimal number of servers to minimize the sum of service costs and delay 
costs? These questions are of much interest to employers and thus to the modeler as well. 
We will explore these questions along with the relevant theory of the M/M/s model (recall 
that s is the number of servers). If 7 customers are present in the system and j < s, then 
all 7 customers are in service and s—j servers are idle. However, if 7 > s, then s customers 
are in service and j — s are waiting in the queue. If a person arrives and a server is idle, 
he enters service immediately; otherwise, he joins the queue of customers awaiting service. 
To model this system as a birth-death process, we remember that the birth rate in state j 
is the arrival rate A; i.e. 4; = A for 7 > 0. If 7 servers are busy, service completions occur 
at arate of ju. Thus, if 7 customers are present, min (j, s) servers will be busy, and so, 
fj = min (j,s) pw. Therefore, the 14/M/s model may be described as a birth-death process 


with parameters 


Aj = A forj>0 
bj = jp for 7 =0,1,...,8 
fj = sp forj>s. (1.2) 
Define 
_A 
Oe Op 


as the traffic intensity. To ensure the stability of the system, p must be less than one. 
That is, the arrival rate must be less than the maximum service rate, which is achieved 
when all s servers are busy. Otherwise, the queue will tend to grow indefinitely in time 
and the system will “blow up” [WA04]. 

Then for p < 1, we substitute (1.2) into (1.1) and obtain the following steady-state 
probabilities: 


For j = 1,2,...,5, we have 


For 7 > s, we have 


se 


Tj 


ToCj 

ToAoA1 Scene Aj—1 
Mip2+ ++ by 
To gf 

3! (si-SpJ) 3f 

(sp) mo 

s!(sI-8)" 


(1.3) 


Note that the steady-state probabilities must add to 1 since at any given time, the system 


is in some state. Thus, 


and so 


e 
ll 


CoO 
yr 
j=0 
s-—1 (oe) 
2 nee 2 i 


She To + caer 3 


co 
Son; =1, 


j=0 


spy 88 A, 
ay 6) or 
j=s 


s—1 


a, 


(oo sfl1 wa, 
—- >" 
1l—p d 


which implies 


Therefore, 


j=0 p 
s—1 j Ss 
8 s 
—? eer 5. (ee) 
j=0 J: 5: ( i, p) 


a ae (1.4) 
& ri + si(l~p 


From equation (1.3), we see that the steady-state probability that all servers are busy (i.e. 


j > 8) is given by 


Now denote the expected number of customers waiting in the queue as Lg. 


P(j 28) 


ll 
“(3 
aad 


= ee (1.5) 


If 7 customers 


are present and j < s, none of the customers need to wait in the queue. However, if j > s, 


there will be 7 — s customers waiting in the queue. Thus, 


ly = So (j- 8) 75 
jas 
[oe 


Changing indices, we obtain 


Hence, from (1.5), we acquire 
(1.6) 
Little’s Formula 


For-any queueing system in which a steady-state distribution exists, the following 
relations hold [Win94]: 


L = \W 
DL, = Ws. 


Note that L = Lg+Ls, and W =W,+W,. Substituting (1.6) into Lg = \W, and solving 
for Wg: 


m = 2 
_ Pies) 
su. (1 — p) 
_ Pes) 
= a (1.7) 


Determining W, and L, along with applying Little’s Formula, we may obtain L and W as 
follows: Since W, = 1/yp, we get 


Thus, 


10 


Then, since L = AW, 
_PG>s) 1 


su-A 
Example 1.1. An employer wishes to determine the number of servers that should work 
on Tuesdays. He approximates a delay cost of $0.05 for every minute a customer waits in 
line. On the average, three customers arrive per minute and it takes a server two minutes 
to complete the service. Each server is paid $12.00 per hour. Interarrival times and 
service times were found to be exponential. How many servers should the employer have 


working on Tuesdays to minimize the sum of service costs and delay costs? 


We have 


’ = 3 customers per minute, and 


LB 0.5 customers per minute. 


Since p = sy must be less than one, 


implies 
s > 6 servers. 
Thus, the employer must have at least seven servers to ensure the stability of the system. 


Then, for s = 7,8,... we will compute: 


Total expected cost _ Expected service cost | Expected delay cost 
Minute a Minute Minute 


Note that each server is paid 12,/60 = $0.20 per minute. Hence, 


Expected service cost 


Minute = $0.20s. 
Also, 
Expected delay cost Expected customers Expected delay cost 
Minute 7 Minute Customer 
= 3(0.05W,) 


0.15W,. 
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Then for s = 7, p © 0.86, from (1.4) we get 7 + 0.0015. Thus, using (1.5), P(j > 7) = 
0.62. From (1.7), Wg = 1.24 minutes. Hence, for s = 7, 


Expected delay cost 


~ $0.19 
Minute $0.19, 


and so 
Total expected cost 


Minute 


Computing the service cost per minute for s = 8 servers in a similar manner, we obtain 


os $1.59. 


$1.60. Therefore, the total expected cost per minute for eight servers cannot possibly be 
lower than that of seven servers since the service cost alone for eight servers is more than 


the total cost for seven servers. Thus, having seven servers is optimal. 
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Chapter 2 


The Erlang Distribution 


When the interarrival times do not fit an exponential distribution function, they 
can often be modeled by an Erlang distribution [Win94]. Let T be a continuous random 


variable with Erlang density function 


R(Rt)k-1e-Rt 
i) - 


where R > 0 is the rate parameter, k is the shape parameter (a positive integer), and t > 0. 


The mean and variance are 


k 
E(T)= R 
and 
k 
var(T) = on 


respectively [Win94]. If k = 1, then f(t) = Re~* (for ¢ > 0), which is an exponential 
distribution with parameter R. Thus, the exponential distribution is a special case of the 
Erlang distribution. More specifically, it is an Erlang type 1. If the interarrival times do 
not seem to fit an exponential distribution, we often consider an Erlang distribution with 
rate parameter kA, shape parameter k, and mean 1//. This can provide greater flexibility 
by being better able to represent the real world [GH85]. 

The sum of k independent, identically distributed exponential random variables 
with mean 1/kyp is an Erlang type k distribution [Win94]._ The Erlang distribution can 
be used to describe queueing models where the service may be a series of k identical phases 


but does not limit its applicability to situations where there are actually phases of service. 
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Note that all phases are independent and identical. Also, only one customer at 
a time can be in the service process. That is, a customer enters service in phase k, then 
phase k—1, and soon. After phase 1 is completed, the customer leaves the system. A 
customer must complete all phases of service before the next one may enter phase k of 


service. 


2.1 The M/E;,/1 Model 


Recall that M/H;,/1 stands for a queueing system with exponential interarrival 
times and Erlang-k service times, where k is the number of phases of service. Even if the 
queueing system does not have phases of service, it is convenient to analyze it in this way 
since each phase may be considered an exponential random variable which allows us to use 
the Markovian property. 

Let pni(t) be the probability of nm customers in the queueing system with the 
customer in service being in phase i, where i = 1,2,...,k. The first phase of service is 
phase k, the second is phase k — 1, and so on; the last phase of service is phase 1. A 
customer leaving phase 1 leaves the system all together. | 


We can write the following set of difference equations: 


Pni(t+ At) = pna(t)(1—AAt— kpAt) + prspa(e)kpAt 
+pn-1i(t)AAt (n>2;1<i<k—1) 
 pnye(t+ At) = pyy(t)(l — AA‘ — kpAt) + pn4ir(t)kyaAt 


+Pn—1,K(t)AAt (nr 2 2) 


pii(t+ At) = pig(t)(l — AAt— kpAt) 
t+priipi(t)kpAt (n=1;1<i<k-1) 


Pik(t+ At) = prz(t)(l— AAt — kpAt) + poi (t)kpAt 
+po(t)AAt 


pott + At) = po(t)(1 - AAt) + P11 (t)kuAt, 


7 whens the probability that an arrival (service completion of phase j) occurs between tine 
—tand time t + Ati is equal to AAt + o(At). (kpAt + o(At)). Note that Atis is an incremental 
element and o(At) becomes negligible when compared to At as At — 0, ie. 


m 24) o(At) _ 
At—0 At 


=0.. 3 ie arene i 
Therefore, all o(At) terms. were ignored i in , the above diffrence equations (GH85]. 
—Let-us first consider the case when n> 1. There are several ways for the system 
to get to state n, j &n customers in | system with the customer in service being i in phasg 4, 


~n>landj=1,2,. a: The system might have been in state n,j at t and had no net 


change during At (nat is, no arrivals and no departures of phase 4). Or the system may : 


| have found itself i in state (n - ay; 3 and had an arrival, or in state n,(j+1) and had a 


. service completion of phase j ie is Tf j= a - k, another possibility i is that the system was in 


state (n +1),1 and had a departure of phase 1 (and thus exiting the system). 
The difference equation for n = 0 may be interpreted as follows: We consider 
~ how the system. may get to state 0 ibe. customers in the system) at time t + At. The 


system might have been i in state 0 at t and had no arrivals during At, or the system might ® 
: have found itself in state 1 and had a service completion (therefore exiting the system), 


~ The corresponding differential-difference equations are found by taking Pn a(t) to 


the left hand side, dividing through by: At, and taking the limit as At — 0. Using the 


: ‘definition of the. derivative, we obtain 


dpni(t) 

Boel a a6, % Ku)Pn i(t) + FuPn, ald: + Nona a(t) 
(n> B1Si< k= 1). 

dpnpe(t) - 

os) s (A+ FiPna(t) + Bets a(0) + Pa at) | 
GS 2 
d alt ‘ : 
spilt) = —(-+ hues) + Rapa sea | | 
(i S$ 4Sh-2) ee cee. Esha, | : 

dpik(t) eg | 

ox as ~(0- helt) + an) + 2000 eI An eg Ge 
er t 
pole = =dpo(t)+ kupra(t). | 


We want to find the steady-atiste difference équations: Thus, we set dpn,«(t)/dt = 0, and 
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so the steady-state difference equations are 


0 = —-A+kp)pni + kupniti t+ rApn-14 (n> 2;1<i<k—-1) 

0 = —(A+kp) pak + kupnsija +APn—1,k (n > 2) 

0 = -(A+kp) pis thupiig: (l<i<k-1) 

0 = —(A+ku)pi~ + kup + Apo 

0 = —Apo + kypij- (2.1) 


Let the random variable T, represent the time spent in the queue, where E({T,] = 
W,. We remember that the expected time to complete each of k phases of service is 1/ ku 
and thus the total expected time for a full service completion is 1/. Thus, the expected 
time a customer spends waiting in line is equal to the expected time for Ng customers in 


line plus the remaining service time of the customer in phase I. 
We = Blt) = BIN + BU, 


1 2 s 
= Piit Piz: + Pik 
ku kp 


kp 

k+1 k+2 k+k 
ce we kis P21 + iF = a ee a ——P2,k 
2k+1 2k +2 2k+k 
+E Patt a Pa tt Paik 
+::++p0 

- ea k(1—1) +2 k(1—1) +k 

= ku Pit ku Piatt kp Pi,k 
k(2—1)+1 k(2—1) +2 k(2—1) +k 
+ fa 2,1 + re p2,2 + + Ris P2,k 
k(3—1)+1 k(3—1)+2 k(3—1) +k 
ar kis D3,1 + ka Bay 8 ee ee 
eRe 

k(n—1) +i 

= Sy [Pee ~ | Pas +o (2.2) 

n=1i=1 


where Ng is the random number for the number of customers in the queue and J is the 
random number representing the phase of service in which the customer being served is 
in. . 

We wish to derive the expected time a customer waits in the queue W, using (2.1). 


Multiplying the first equation in (2.1) by z*”—)+4, the second by z*", the third by z*, the 
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fourth by z*, and the fifth by z° gives 


0 = -(A+ ky)zh—D+4tp,, + kyzklr tin, iy + AzkO—Dtty, 4 where 
(n > %1<i<k-1) 
= —(A+kp)z* pap + kuz’ pags + AZ" Paik (n= 2) 
—(A+ ku)z'piat ku2'riigt (L<i<k-1) 
= —(A+kp)z*pip + kuz*po1 + Az" po 


ooo Oo 
Il 


= —2po + ky2°pi1. (2.3) 
Let 


Gz) = zpiatz*pe+...+2* pip 
taht + zh +... +27 poe +... + Do 


co k 
AO Hons + Po- 


n=li=1 


Expanding each equation in (2.3) and summing up all terms, we have: 


0= : (G(z) — po) - (1 + =) G(z) + po + ate). 
Let r=A/kp. Then, 
0= - (Ga) =p) = C49 OG) Fo PSG): 


It follows that 


= po(1 — 2) 
CU) qa: a r) + rzktl? 
where 
po=i- - 
Thus, 
@iye= [t—2(ltr)+rz*1] -(1-z)[-(l+r) + (k4+1)r2*] (1- *) 
[l—2(l+r)+rz¢41/? ue) 


To find G’(1), we have to use L’Hopital’s rule twice. We obtain 
. d 
he 


G'(1) = 
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From (2.2), we conclude 


Since the total wait time is a sum of the expected wait time in queue and the expected 


wait time in service, we have: 


Ww = W,+W, 
1 
= sd Ter 
(k+1)\+2k(u—2) 
2ku (uw — A) 


Using this together with Little’s Formula, we also obtain 


_ (k+1)»# 
ky (ps — 2) 
and 
Lb = \w 


(k+1)\_ 1). 
‘ (sere si 7 
— y(RAAAt 2k (uA) 
= a( 2kp (wu — A) ). 


2.2 An Example for M/E;/1 


We wish to determine whether the marketing department should rent a slow or 
fast copy machine. The department estimates each employee’s time to be worth $15.00 
per hour. The slow copier rents for $4.00 per hour, and it takes an employee an average 
of 10 minutes with a variance of 19 minutes to complete copying. The faster copier rents 
for $8.00 per hour, and it takes an employee an average of 6 minutes with a variance 


of 9 minutes to finish copying. An average of four employees per hour need to use the 
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machine. The interarrival times were found to be exponential. Which machine should 
the department rent to minimize the total expected cost per hour? 


Let us first analyze the slow copier. Note that 


A = 4 employees per hour 


and 
1 : 
— = 10 minutes per employee 
1 
mane hours per employee, 
so 
1 ; 
js = — employees per minute 


10 
= 6 employees per hour. 


Let us check that the copy times satisfy 
var T < [E(T)]’. 


Since 19 < 10?, an Erlang distribution can be fitted to the copy times [Win94]. Let us 


determine k. 


“ai 
var T = ge 
1 
; k (76) 
k = 5.26. 


Thus, k =5o0rk=6. For k = 5, var T = 20, and for k = 6, var T = 16.67. Hence, 
we choose k = 5 since the observed variance is closest to that achieved with k = 5. 
Therefore, the copying times fit an Erlang distribution with rate parameter ku = $ and 
shape parameter k = 5. We have an M/Es5/1 model. We must also check system stability 
[Win94], ive. 
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thus the queueing system is stable. 
Now we carry out the calculations (using hours as the time unit) to determine the 


total expected cost per hour for the slow copy machine as follows: 


Total Expected Cost _ Expected Service Cost rs Expected Delay Cost 


Hour Hour Hour ; 
where 
Expected Delay Cost _ (Ss | ( Delay ) 
Hour Hour Employee , 
Thus, 


Total Expected Cost 


=4+4(15W,). 
Hour TSS) 


Note that the expected delay cost per hour is the employee’s wage per hour. Calculate 


W, as follows: 


(k+1)A 
a 2kp (wu — 2d) 
_ $6) 
~  60(2) 
1 
=e hours. 
Hence, for the slow machine 
Total Expected Cost 60 
i = sh 
Hour Bs) 
= $16.00. 


Thus, the slow copier has an expected cost of $16.00 per hour. 


Similarly, we may analyze the fast copy machine. We have 


A = 4 employees per hour 
1 
ji = 6 minutes per employee 
1 
= 59 hours per employee 
iL 
fs = = employees per minute 


6 
= 10 employees per hour. 
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Since var T < [E(T)]’, or 9 < 6”, an Erlang distribution may be fitted to the copy times. 
We also verify that p = A\/p = 4/10 < 1, so the system is stable. To find k: 


var T = a 
1 
9 = 
1\2 
k (3) 
k= 4 
Similar to the case for the slow copy machine, 
Total Expected Cost _ Expected Service Cost re Expected Delay Cost 
Hour = Hour Hour 
= 8+4(15W,). 
Calculate W, as follows: 
(K+1)A 
W, pa ec ies Et 
7 2ku (ud) 
_ 5) 
~~ 80(6) 
= = hours 
ne hie 
Thus, for the fast copy machine, 
Total Expected Cost __ 8-4 60 
Hour 7 24 
= $10.50. 


Therefore, the fast copier has an expected cost of $10.50 per hour. 

Since the total expected cost per hour is $16.00 for the slow machine versus $10.50 
for the fast copier, the marketing department should rent the fast copy machine to minimize 
the total cost per hour. 

Note that we may also determine the expected number of customers in queue (Lg) 
for the slow and fast copy machines using the steady-state formula that we derived. The 
slow copier gives Lg = 4/5 employees, and the fast copier gives Ly = 1/6 employees. Of 


course, the smaller expected wait time produced a smaller expected number in queue. 
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Chapter 3 


Model Selection 


Finding the appropriate probability distribution that describes a set of data is 
imperative for a model to produce meaningful results. For the queueing systems that we 
are interested in, determining whether the actual data are consistent with the assumption of 
exponential interarrival times and service times is essential. Goodness of fit tests are often 
used to test a set of data for fitting a probability distribution. The most common tests are 
the chi-square goodness of fit test, Kolmogorov-Smirnov test, and Anderson-Darling test. 
In this chapter, we will demonstrate how to use these three tests to choose appropriate 
probability distributions for the given data set. 

All of these tests compare the data to a theoretical distribution (with the popu- 
lation mean m replaced by sample mean ™). Given 7m, we wish to conduct a hypothesis 
test to determine whether 1}, to, ...,t, represent a random sample from a random variable 


with a given density function f(t). That is, we want to test the following hypotheses: 


Hy : t1,te,...,t, is a random sample from a random 
variable with a given density function f(t), 
Ha : t1,t9,...,¢, is not a random sample from a random 


variable with a given density function f(t). 


Let us consider queueing models in which we want to fit to an exponential distribution. 
Suppose the interarrival times of a queueing system have been observed to be #y, ta, ..., tn, 


where ¢; is the time between the (i—1)st and ith arrival. We want to determine an estimate 


22 


of the arrival rate \ from the observed data. A common method is using the maximum 
likelihood estimate [KPW04]. 


Definition 3.1. The likelihood function is 


n 
L(6) = | [ PCT; = ty | 9) 
j=l 
and the maximum likelihood estimate of 0 is the value that maximizes the likelihood 
function. 


This estimate is found by setting the objective function L(@) and then determin- 
ing the parameter value that optimizes the function. We will show that the maximum 


likelihood estimate of the arrival rate is given by 


n 


A = nn 
dit 
i=1 
and the mean of interarrival times 1A = 6 can be estimated by 
n 
dit 
=F = 
n 


By the definition of the likelihood function, 


n 
L6) = [[PG=4 1) 

j=l 

= [fe 14) 
j=l 

n a 

= [[e te # 
j=l 

—y 

= 9%e =? 


Note.that when calculating the likelihood function for a single point, we interpret 


P(T; =t; | 0) as f(t; | 0); ie. we use the density function. 
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Instead of maximizing L (@), we will maximize / (6) = In L (6). This can be done 


since both L (6) and 1 (@) have their maximum point at the same value. Then, 


nr 
St 
l(@) = —nind — = 
3 
tj 
Wp) — 2 1 _ 
0) = -7+4-=0 
nr 
0 = —nd+ 5 >t 
i=1 
n 
dit 
621 =F 
n 2 


3.1 Chi-Square Goodness of Fit Test 


The first step in a chi-square goodness of fit test is to break up the set of in- 
terarrival times into k adjacent intervals [ao, a1), [@1,@2),.-.,[@e-1,@%), where a, = oo. 
Assuming that f(t) governs interarrival times, we determine the number of expected in- 
terarrival times, denoted by ej, that are in interval 7. To do this, compute the expected 
probability p; of the data that would fall in the jth interval, where 

p= >, Ba) 
3-1 tj Sa; 
and p is the mass function of the fitted distribution. Hence e; = np; is the expected 
number of the data that will fall in the jth interval. Denote the number of observed data 


in the jth interval [a;_1,a;) by 0;. The chi-square test statistic is 


2 _ wr (oy - 63)? 
v= y oa 
A small test statistic value demonstrates a good fit. This occurs when the observed 
interarrival times are near the expected interarrival times. 

Given a value of the desired Type I error (rejecting H, when it is true) a, the 
critical value is Nee where 7 is the number of parameters that must be estimated 
for the interarrival time distribution. Note that r = 1 for the exponential distribution. 


Thus, accept H, if x? < Xk-2,1-a" 
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Accuracy of a test i8 an important factor in modeling. The requirements for 
accuracy of the chi-square test are not concrete, but there are a few guidelines that are 
often followed. The intervals should be chosen such that p; = po = ... = px; this is called 
the equiprobable approach. This is satisfied for 3 < k < 30 intervals and e; > 5 for all j 
[Win94], [LK00]. The lack of rigid rules for interval selection for the chi-square goodness 
of fit test is one of its major drawbacks. 

However, for the equiprobable approach for interval selection, the chi-square test 
is said to be unbiased since it is more likely to reject a false null hypothesis than a true one. 
If many e;’s are small and not equal, a highly biased yet valid test is possible. According 
to [LK00], the equiprobable approach along with the conditions that e; > 5 for all j and 


k > 3 guarantees a valid and unbiased test. 


Example 3.1. The interarrival times (in minutes) observed at a bank are as follows (in 
ascending order): 0.1, 0.2, 0.2, 0.8, 0.3, 0.5, 0.6, 0.6, 0.7, 0.8, 0.9, 1.2, 1.5, 1.6, 1.8, 
2.2, 2.8, 8.0, 3.5, 3.8, 3.9, 5.3, 6.1, 9.4, 11.0. Use the chi-square goodness of fit test 
with a = 0.05 to determine if it is reasonable to conclude that the observations follow an 


exponential distribution. 


There are 25 observations with a mean arrival rate 


A= a = 0.405 arrivals per minute. 
dot 
i=1 


Now test whether the data are consistent with an exponential random variable, say A, with 
density f(t) = 0.405e~°-495, Choose k = 5 intervals to ensure that the probability that an 
observation from A falls into each of the five categories is pj = 0.20 forl <j <5. Thus 
e; = 25(0.20) = 5 for each interval. Then set the interval boundaries. We must first 


calculate the cumulative distribution function, F(t), for A: 


F(t) = P(A<t) 


t 

‘ 0.405¢7 0:4052 gap 
0 

eee, ea e 0-405t 


Then the five intervals will be [0, a1), [a1, a2), [a2, a3), [@3, 24), [a4, 00), where 


F(a1) => P(A < ay) = 0.20 


25 


F (a2) = P(A < ag) = 0.40 
F(a3) = P(A < a3) = 0.60 
F (a4) = P(A < a4) = 0.80. 


From F(t) = 1 — e~ 9-4, we can use natural logarithms to solve for t as follows: 


F(t) = p= e0-405t 
Ine~°-45% = In(1— F(t)) 
a In(1 — F(t)) 
a —0.405 
In(1 — F(a;)) : 
a; 405 for all 7. 


It follows that 
a, = 0.55, ag = 1.26, a3 = 2.27, and a4 = 3.98. 


The number of observations in each interval is 
01 = 6, 09 = 6, 03 = 4, 04 = 5, and os = 4. 


Thus, 


_ (6-5)? (6-5)? (4-5)? (5-5)? (4-5)? 
ee Sih ee, Se ee er 


= 0.20 + 0.20 + 0.20 + 0 + 0.20 
= 0.80. 


The critical value is X3.0.95 = 7.81, and since x? = 0.80 < 7.81, we cannot reject the 
null hypothesis. Hence, with 95% certainty, the null hypothesis is not rejected and the 
exponential distribution with 1 = 0.405 arrivals per minute is a plausible model for the 


interarrival times. | 


3.2. The Kolmogorov-Smirnov Test 


The Kolmogorov-Smirnov (K-S) test is another method that can be used to assess 
whether the observations T},7>,...,T;, are an independent sample from the exponential 


distribution. The K-S test compares an empirical distribution function with the fitted 
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distribution function. This test has advantages and disadvantages. ‘To perform the test, 
the data does not need to be grouped, and hence fhe problem of interval specification is 
eliminated. Although the original K-S test requires all parameters to be known (i.e. the 
parameters should not be estimated from the sample data), applying it for discrete or for 
any continuous distribution with estimated parameters yields a conservative test. That is, 
the actual probability a/ of rejecting the null hypothesis when it is true is at least as small 
as the stated probability a [LK00]. 

For the K-S statistic, define the empirical distribution function F,(t) as the right- 
continuous step function where Fy, (Tuy) =i/n fori =1,2,...,n. Also, the distribution 
function F(t) is assumed to be continuous over the range of data. A natural assessment of 
goodness of fit is a measure of closeness of the empirical and fitted distribution functions. 
In the case when all parameters are known, the test statistic is the largest vertical distance, 
Dn, of F,(t) and F(t) for all t, where Dp, is defined by [KPW04] 


Dr = sup |Fn(t) — F(#)I. 


However, since in our case not all parameters are known (we estimate the mean from the 
observations), we will use the adjusted test statistic 
0.2 0.5 
Dn - — 0.26 + —}. 
Pe) Aes a) 
Once more, we do not reject H, if this test statistic is small. A commonly used critical 


value for this situation is c}_, = 1.094 for a = 0.05 [LK00]}. 


Example 3.2. Use the K-S test with a = 0.05 on the previous example to determine if it 


is reasonable to conclude that the observations follow an exponential distribution. 


We first set up a table to compute F,,(t), F(t), and |F,(t) — F(¢)| foreacht. Then 
we determine D,, from the table (see Appendix A). Since the critical value is cp.95 = 1.094 


and 


0.2 0.5 
Dos — — | | V25 + 0.26 + —=] #0. < 1.094 
( 25 a) ( + +) 0.6798 < 1.094, 


again we cannot reject the null hypothesis. Hence, using the K-S test with 95% certainty, 
the null hypothesis is not rejected and the exponential distribution is a plausible model for 


the interarrival times. O 
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3.3. The Anderson-Darling Test 


The Anderson-Darling (A-D) test is a similar test to the K-S test, however it 
detects discrepancies in the tails. This is an important characteristic since many distribu- 
tions mostly differ in their tails [LK00]. Denote the A-D test statistic for the case when 


all parameters are known by 
[o.e) 
Aan [ F(t) - FOP vOFE at 
00 


where the weight function 
1 


U() = eae 
= FOE -FO 
Then A? is the weighted average of the squared differences between the empirical and 
model distribution functions. The weights are largest for F(t) close to 0 and 1 (the left 
and right tails, respectively). Let Z; = F(Ty)) for i =1,2,...,n. Then the test statistic 


for carrying out actual computations can be written as 


n 
-_ {> (2% = 1) [In Zit In(1 = Zni-3)]} 
AZ = i=1 
ig n 
The values for F;,(t) and F(t) are calculated the same way as for the K-S test. 


—-n. 


Since A? is a weighted distance, we want to reject the null hypothesis if A? is too 
large. Once again, an adjusted test statistic is available that takes into account that the 
mean was estimated from the data. Thus, reject H, if 

(1 + =) A? > 1.326 
for a = 0.05 [LKO0}. \ 
Example 3.3. Use the A-D test with a = 0.05 on the previous example to determine if it 
is reasonable to conclude that the observations follow an exponential distribution. 

We set up a table to compute F,,(t), Z;, and 241-4 for each t (see Appendix B). 
Then compute A2 s 0.4587. Since the critical value is do.95 = 1.326 and 


(1 + ss) Ad, ~ 0.4697 < 1.326, 


again we cannot reject the null hypothesis. Thus, using the A-D test with 95% certainty, 
the null hypothesis is not rejected and the exponential distribution is a plausible model for 


the interarrival times. oO 
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3.4 Concluding Remarks 


We first note that any model is only an approximation of reality, however the 
model may still be useful. While selecting a distribution model, it is very important to 
follow the principle of parsimony. The principle states that unless there is strong evidence 
to use a more complex model, it is preferred to choose a simpler model. Most hypothesis 
tests use this principle as well. Unless there is strong evidence to do so, do not reject 
the null hypothesis (and thus claim a more complex model of the population needs to be 
found). Finally, it is wise to keep focused on the problem that is to be solved rather than 
spend abundant energy searching for the perfect model. The reasoning behind this is that 
if a complex model is very close to the observations, there is no guarantee that the model 
matches the population from which the data were sampled [KP W04]. 

There are two main approaches to model selection; they are the judgment-based 
and score-based approaches [KPW04]. The judgment-based approach allows for a decision 
to be based on the success of particular models in similar situations or on how well the 
data compare to the empirical distribution using a graph. 

The score-based approach assigns each model a score and the model with the best 
score is chosen. Common scores include the lowest value of the Kolmogorov-Smirnov test 
statistic, the lowest value of the Anderson-Darling test statistic, the lowest value of the 
chi-square goodness of fit test statistic, the highest p—value for the chi-square test, and the 
highest value of the likelihood function at its maximum. Overall, the analyst’s judgment 
is required at the very least in deciding which algorithm to choose. 

Each hypothesis test has its advantages and disadvantages. The chi-square good- 
ness of fit test is commonly used since it may be applied to any hypothesized distribution. 
Also, the critical value of the test is easily adjusted depending on the number of param- 
eters estimated from the data. A valid chi-square test using the equiprobable approach 
for interval selection always produces an unbiased test. A test is said to be unbiased if it 
is more likely to reject a false null hypothesis than when it is true. However, the major 
drawback of this test is the lack of rigid rules for interval selection. In some cases, different 
choices of intervals can lead to different conclusions [LKO0]. 

For the Kolmogorov-Smirnov test, the data does not need to be grouped, which 


eliminates the problem of interval specification. However, this test requires an adjusted 
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test statistic when parameters are estimated from the data. The critical values are not 
readily available for discrete data and must be computed from difficult formulas. Lastly, 
the K-S test statistic gives the same weight to the difference |F;,(t) — F°(t)| for all values of 
t, but many distributions differ primarily in their tails [LK00]. 

The Anderson-Darling test has an advantage since its test statistic is a weighted 
average where the weights are largest in the tails, unlike the K-S test statistic. Similar 
to the K-S test, the A-D test statistic must be adjusted for the case when parameters are 
estimated from the data. The critical values for this test also must be computed from 
difficult formulas for discrete data [LK00]. 
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Chapter 4 


Simulation 


Simulation is a common method of analyzing queueing systems as an alternative 
method to the analytical approach. Perhaps the biggest advantage of simulation is its flex- 
ibility in the sense that it is possible to create a program for an individual queueing system 
without needing as many simplifying assumptions as the analytical method. Queueing 
systems with non-exponential service times, a limited capacity waiting room, and many 
other situations which are extremely difficult to analyze directly can be explored using 
simulation. Simulation should be used when an analytical solution may not be found or 
it does not have acceptable approximations [GH85]. 

Another advantage of simulation is that the user sees the action through time. 
Creating simulation models to see the effects of changing system parameters is much more 
cost efficient: than to change the system in real life. There are three basic phases of 
simulation: data generation, bookkeeping, and output analysis. In queueing systems, 
the first phase generates the interarrival and service times with appropriately selected 
probability distributions. The biggest problem in queueing simulation is bookkeeping. 
Bookkeeping involves keeping track of arrivals, departures, busy and idle servers, queue 
length, clock time, and status of the server. In other words, timing and bookkeeping is a 
challenge [WA04]. The outputs are then analyzed to determine which parameters of the 
queueing system should be changed to optimize system performance. 

The model that we will consider is a single server queueing system. Customers 
arrive at random times to the system. If a customer arrives and the server is busy, the 


customer joins the end of a single queue. The system starts empty and idle. We then 
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simulate the system for a user-defined length of time, called the close time. At this time, 
no new arrivals are allowed to enter the system, but customers already present are served. 
The simulation terminates when the last customer leaves. This model assumes that the 
times between arrivals are exponentially distributed and the service times are Erlang type 
k. The purpose of the simulation is to collect statistics on the system behavior in order 


to optimize its performance. 


4.1 Description of the Simulation 


Visual Basic for Applications (VBA) was used to simulate the queueing system 
discussed in this chapter. The VBA code takes care of all the timing and statistical 
bookkeeping as the simulation runs [Alb00]. The program code may be found in Appendix 
C. The key idea of the simulation is one of scheduling events. At any given time, there is 
a list of scheduled events of two types. The first type is an arrival. Each time a customer 
arrives to the system, the next arrival is scheduled at some random time in the future. The 
second type of event is a service completion (or departure). Each time a customer goes 
into service (possibly after waiting in the queue), a departure is scheduled for a random 


time in the future. 


Generating Exponential and Erlang Type k Service Times 


We first want to generate arrival times from the exponential distribution with 


parameter . We recall that the probability distribution is given by 
F(z) =1-—e"** (¢>0). 


We use VBA’s random number generator to produce a uniform-(0,1) random number r. 
Solving the equation 


r=1—e* 


for t is referred to as the analytical inversion process of generating random variables from 


an exponential distribution [GH85]. The above equation becomes 


l-r=e 
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Note that it does not matter whether we use 1—r or r since both are uniform-(0,1) random 


numbers. Thus, we may write 


and by taking natural logarithms of both sides, we obtain 


lnr 


t=—-— 


a 
This is the formula used in our VBA program to generate random exponential interarrival 
times. 
To generate Erlang random variables, we may use the fact that an Erlang type k 
random variable, z, is the sum of k independent, identically distributed exponential random 
variables with mean 1/ky [Win94]. Thus, using VBA’s random number generator to 


produce uniform-(0,1) random numbers rj, r2,...,7~, we have 


8 

ll 
Me 

| 
=| 
2 {5 
Ne, 


which is the formula used to generate random Erlang type k service times. 
We will now explain the major steps of the program after the user enters the 
following inputs: the customer arrival rate, the mean service time per customer, the number 


of phases of service, and the closing time. 


Step 1: 
The program begins by initializing the system. Set the clock time to zero, set the status 


of the server to idle, and schedule the first arrival. Go to step 2. 


Step 2: 

Determine whether the next event will be an arrival or a departure as follows. Find the 
minimum of the scheduled event times. If the next event time is that of an arrival, reset 
the current clock time to the time of the arrival and go to step 3. If it is a departure, 
reset the current clock time to the time of the departure and go to step 4. In either case, 
increase the wait times of everyone in the queue by the elapsed time between the previous 


and next event. 
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Step 3: 

Check to see if the server is busy. If yes, put this arrival at the end of the queue, keeping 
track of the arrival time (for later statistics). If the server is idle, place this customer into 
service and schedule his departure. Schedule the time of the next arrival. If this time is 
after the closing time, do not allow the next arrival and do not schedule any future arrivals. 


Go to step 5. 


Step 4: 

If there is at least one customer in the queue, send the customer from the front of the queue 
to the server and record his wait time for later statistics. Move all other customers up one 
space in the queue and schedule a departure. If there is no queue, set the server’s status to 
idle; do not schedule a departure event. Increase the number of served customers by one. 
Go to step 5. 


Step 5: 
If the clock time is greater than the close time and the server is idle, calculate the outputs 
and terminate the program. Otherwise, go to step 2. 


The flow chart may be found in Appendix C. 


4.2 The M/E;,/1 Example Revisited 


Recall that in Chapter 2 we considered an example where we used analytical 
methods to determine summary measures for the M/F;/1 model. Here, we want to 
compare the theoretical values of the expected wait time in queue (W,) and the expected 
number of customers in queue (Zq) found in Example 2.2 with the corresponding values 
using a simulation model. 


From Example 2.2, we have for the slow copy machine 


= 4 employees per hour, 


i 
= 's hours per employee, and 


zwei y 


= 5 phases of service. 


er 3 
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Theoretical values for the expected wait time and the expected number in queue were 


W, = = hours 


Lg : employees. 


For the fast machine we have 
bas ae hours per employee 
k = 4 phases of service 


1 
Wy = 94 hours 
Ig = ; employees. 


We use the VBA program to simulate the use of each copy machine. 


Inputs and Outputs 


There are four inputs for each copier: the employee arrival rate 4, the mean service 
time per employee 1’, the number of phases of service k, and the closing time. At the 
end of the simulation, we want to display the summary measures, which include: average 
time and maximum time in queue, average number and maximum number of employees in 
the queue, and the fraction of time that the server is busy. Also, the simulation program 
outputs the number of employees processed and the probability distribution of number in 


queue. 


From [Ros02}, we determined to perform one hundred runs for the slow and fast 
machines. Each run simulates an eight hour period. Also, by [Win94], we expect the 


average of the simulation values to be close to the steady-state values found in Chapter 2. 


Slow Copy Machine 


We calculated the average wait time and the number in queue to be 0.20055 hours 
and 0.80102 employees, respectively. Both averages are very close to the steady-state values 


of 0.20 hours and 0.80 employees, respectively. 
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Fast Copy Machine 


We calculated the average wait time and the number in queue to be 0.04158 
hours and 0.16644 employees, respectively. Again, both averages are very close to the 


steady-state values of 0.041667 hours and 0.166667 employees, respectively. 


In summary, simulation allows the user to obtain a realistic analysis of the queue- 
ing system. This method is commonly used when a queueing model is extremely difficult 
or impossible to analyze analytically due to the complexity of the queueing system. Simu- 
lations allows us to create hypothetical situations without costly real-time experimentation 


nor possible loss of clientele due to long lines and long wait times. 


Ng 


Appendix A 


Kolmogorov-Smirnov Test 


Example 
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Table for Kolmogorov-Smirnov Test Example 3.2 


ne 25 
i t 
1 0.1 
2 0.2 
3 0.2 
4 0.3 
5 0.3 
6 0.5 
7 0.6 
8 0.6 
9 0.7 
10 0.8 
11 0.9 
12 1.2 
13 1.5 
14 1.6 
15 1.8 
16 - 2.2 
17 2.3 
18 3 
19 3.5 
20 3.8 
21 3.9 
22 5.3 
23 6.1 
24 9.4. 
25 no 
Mean t = 2.472 
D, = 0.134837225 


Adjusted Statistic = 0.679847524 


Fn(t) =i/n 
0.04 


F(t) 
0.039645771 
0.077719756 
0.077719756 
0.114284267 
0.114284267 
0.183121878 
0.215507641 
0.215507641 
0.246609446 
0.276478195 
0.305162775 
0.384571738 
0.454905506 
0.476516198 
0.517201231 
0.589330957 
0.605612248 
0.702871993 
0.757282632 
0.785021408 

0.7935444 
0.882816353 
0.915215076 
0.977687071 
0.988319543 


| Fr(t) - F(t) | 

0.000354229 
0.002280244 
0.042280244 
0.045715733 
0.085715733 
0.056878 122 
0.064492359 
0.104492359 
0.113390554 
0.123521805 
0.134837225 
0.095428262 
0.065094494 
0.083483802 
0.082798769 
0.050669043 
0.074387752 
0.017128007 
0.002717368 
0.014978592 

0.0464556 
0.002816353 
0.004784924 
0.017687071 
0.011680457 
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Appendix B 


Anderson-Darling Test Example 


Table for Anderson-Darling Test Example 3.3 


n= 25 
i t F,(t) =i/n 2, = F(t) Znati -—« (2-1) [In Z; + In (1-Z,,;.)] 
1 0.1 0.04 0.039646 0.988320 7.677609173 
2 0.2 0.08 0.077720 0.977687 19.07170439 
3 0.2 0.12 0.077720 0.915215 25.1114167 
4 0.3 0.16 0.114284 0.882816 30.19155515 
5 0.3 0.2 0.114284 0.793544 33.72062638 
6 0.5 0.24 0.183122 0.785021 35.58302197 
7 0.6 0.28 0.215508 0.757283 38.35801475 
8 0.6 0.32 0.215508 0.702872 41 .2252672 
9 0.7 0.36 0.246609 0.605612 39.61629166 
10: 0.8 0.4 0.276478 0.589331 41.33622821 
11 0.9 0.44 0.305163 0.517201 40.2163712 
12 1.2 0.48 0.384572 0.476516 36.86610482 
13 1.5 0.52 0.454906 0.454906 34.86154191 
14 1.6 0.56 0.476516 0.384572 33.12064231 
15 1.8 0.6 0.517201 0.305163 29.67862673 
16 . 2.2 0.64 0.589331 0.276478 26.42415051 
17 2.3 0.68 0.605612 0.246609 25.89466682 
18 3 0.72 0.702872 0.215508 20.8354628 
19 35 0.76 0.757283 0.215508 19.26727583 
20 3.8 0.8 0.785021 0.183122 17.32807683 
21 3.9 0.84 0.793544 0.114284 14.45680537 
22 5.3 0.88 0.882816 0.114284 10.57788406 
23 6.1 0.92 0.915215 0.077720 7.627605028 
24 9.4 0.96 0.977687 0.077720 4.86317353 


25 11 1 0.988320 0.039646 2.557911903 


Mean t, = 2.472 


A,” = 0.458721 
Adjusted Statistic 0.469731. 


Appendix C | 


Flow Chart and VBA Code 


Set clock time = 
tlme of next arrival 


YES 


Increase the number 
In queue by 1 


Set time of next arrival 
= Clock time + IT 


Is the server busy? 


Generate exponential 
interarrival time (IT) 


Generate new 
random numbers 


Initlallze counters and 
status Indicators 


YES, process 
an arrival 


NO, process 


Is time of next a departure 


arrival < time of 
next departure? 


Set clock time = time 
of next departure 


Is the number In 
queue = 07 


NO 
Set server status to Set server status Decrease the number 
busy to idle In queue by 1 
Generate Erlang Generate Erlang 
service time (ST) service time (ST) 
Set time of next departure Set time of next departure! 
= Clock time + ST 


Is the clock time > close 
time and the server Idle? 


NO 


YES Print results and 
terminate program 
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Code for the Simulation Program 
Option Explicit 


" Declare system parameters. 

‘MeanIATime - mean interarrival time (reciprocal of arrival rate) 
MeanServeTime - mean service time 
MaxAllowedInQ - maximum number of customers allowed in the queue 
CloseTime - clock time when no future arrivals are accepted 
K - used to generate Erlang number with shape parameter K 
MeanServeTimeDivByK - MeanServeTime/K 
ProdKRand - product of K random numbers between 0 and 1 


Dim MeanJATime As Single, MeanServeTime As Single, _ 
CloseTime As Single, K As Integer, MeanServeTimeDivByK As Single, ProdKRand As Single 


Declare system status indicators. 

NumInQ - number of customers currently in the queue 

ServerBusy - the status of the server, 0 for idle, 1 for busy 

ClockTime - current clock time, where the inital clock time is 0 

TimeOfLastEvent - clock time of previous event 

EventScheduled(i) - True or False, depending on whether an event of type i is scheduled or not, 
for i=0,1, where i=0 corresponds to arrivals and i=1 corresponds to service completions 

TimeOfNextEvent(i) - the scheduled clock time of the next event of type i (only defined when 
EventScheduled(i) is True) 


Dim NumInQ As Integer, ServerBusy As Integer, ClockTime As Single, TimeOfLastEvent As Single, 
7 EventScheduled(Q As Boolean, TimeOfNextEvent() As Single 


‘Declare statistical variables. 
NumServed - number of customers who have completed service so far 
MaxNumiInQ - maximum number who have been in the queue at any point in time so far 
MaxTimelnQ - maximum time any customer has spent in the queue so far 
-TimeOfArrival(i) - arrival time of the customer currently in the i-th place in the queue, for i>=1 
TotalTimeInQ - total customer-time units spent in the queue so far 
TotalTimeBusy - total server-time units spent serving customers so far 
SumOfQTimes - sum of all times in the queue so far, where the sum is over customers who 
who have completed their times in the queue 
QTimeArray(i) - amount of time there have been exactly i customers in the queue, for i>=0 


Dim NumServed As Long, MaxNumInQ As Integer, MaxTimeInQ As Single, _ 
TimeOfArrival() As Single, TotalTimeInQ As Single, TotalTimeBusy As Single, _ 
SumOfQTimes As Single, QTimeArray() As Single, NumLost As Integer, KRand() As Single 


Sub MainO 
Dim NextEventType As Integer, i As Integer, m As Integer 


‘ Always start with new random numbers. 


a err eee 


= 


—_—__—, 


ee 


A 
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Randomize 


' Clear previous results, if any, from the Report sheet. 
“Call ClearOldResults 


' Get inputs from the Report Sheet. 
MeanIATime = 1 / Range("ArriveRate") 
MeanServeTime = Range("MeanServeTime") 
CloseTime = Range("CloseTime") 

K = Range("K") 


‘Compute MeanServeTimeDivByK 
. MeanServeTimeDivByK = MeanServeTime / K 


‘The next two arrays have an element for arrivals (element 0) and one for the server. 
ReDim EventScheduled(2) 
ReDim .TimeOfNextEvent(2) 


"Set counters, status indicators to 0 and schedule first arrival. 
Call Initialize 


' This array has an element for each K; each will be a random number between 0 and 1 
ReDim KRand(K) 
Fori=1ToK 
KRand(i) = Rnd 
Next 


‘Create the product of K random positive numbers less than 1. 
Form=1ToK 
ProdKRand = ProdKRand * KRand(m) 
Next | . 


' Keep simulating until the last customer has left. 
Do 


Find the time and type of the next event, and reset the clock. Capture the index of the finished 


_ "server in case the next event is a service completion. 


Call FindNextEvent(NextEventType) 


‘ Update statistics since the last event. 
Call UpdateStatistics a 


'NextEventType is 1 for an arrival, 2 for a departure. 
If NextEventT ype = 1 Then 
Call Artival 
Else 
- Call Departure 
End If ‘ 
Loop Until (ClockTime > CloseTime And ServerBusy = 0) Or (Not EventScheduled(0) And Not 


EventScheduled(1) And ServerBusy = 0) 
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"Report the results. 
Call Report 
End Sub 


Sub ClearOldResults() 
With Worksheets("Report") 
.Range("B12:B23").ClearContents 
With .Range("A26") 
Range(.Offset(1, 0), .Offset(0, 1).End(xIDown)).ClearContents 
End With 
End With 
End Sub 


Sub [nitializeQ 


"Initialize ProdKRand to 1 
ProdKRand = 1 


‘Initialize system status indicators. 
ClockTime = 0 
ServerBusy = 0 
NumInQ = 0 
TimeOfLastEvent = 0 


‘Initialize statistical variables. 
NumServed = 0 
SumOfQTimes = 0 
MaxTimelInQ = 0 
TotalTimeInQ = 0 
MaxNumInQ = 0 
TotalTimeBusy = 0 


"Redimension the QTimeArray array to have one element (the 0 element, for the amount of time when 
‘there are 0 customers in the queue). 

ReDim QTimeArray(1) 

QTimeArray(0) = 0 


‘ Schedule an arrival from the exponential distribution. 
EventScheduled(0) = True 
TimeOfNextEvent(0) = -MeanIATime * Log(Rnd) 


‘Don't schedule any departures because there are no customers initially in the system. 
EventScheduled(1) = False 


End Sub 


Sub FindNextEvent(NextEventType As Integer) 
Dim i As Integer, NextEventTime As Single 


‘NextEventTime will be the minimum of the scheduled event times. Start by setting it to a large value. 
NextEventTime = 10 * CloseTime 


‘Find type and time of the next (rnost imminent) scheduled event. Note that there is a potential 
‘event scheduled for the next arrival (indexed as 0) and for a server completion. 
For i=0To 1 
If EventScheduled(i) Then 


‘If the current event is the most imminent so far, record it. 
If TimeOfNextEvent(i) < NextEventTime Then 
NextEventTime = TimeOfNextEvent(i) 
If i=0 Then 
' Arrival case 
~ NextEventType = 1 
Else 


" Departure case 
NextEventType = 2 


End If 
End If 
End If 
Next 


‘Update the clock to the time of the next event, if there is one. 
If EventScheduled(O) Or EventScheduled(1) Then 
ClockTime = NextEventTime 
End If 

End Sub 


Sub UpdateStatistics() 
Dim TimeSinceLastEvent As Single 


 TimeSinceLastEvent is the time since the last update. 
TimeSinceLastEvent = ClockTime - TimeOfLastEvent 


' Update statistical variables. 
QTimeArray(NumInQ) = QTimeArray(NumInQ) + TimeSinceLastEvent 
TotalTimeInQ = TotalTimeInQ + NumInQ * TimeSinceLastEvent 
TotalTimeBusy = TotalTimeBusy + ServerBusy * TimeSinceLastEvent 


"Reset TimeOfLastEvent to the current time. 
TimeOfLastEvent = ClockTime 
End Sub 


Sub ArrivalQ 
Dim i As Integer, j As Integer 


" Schedule the next arrival. 
TimeOfNextEvent(0) = ClockTime - MeanIATime * Log(Rnd) 


‘ Cut off the arrival stream if it is past closing time. 
If TimeOfNextEvent(0) > CloseTime Then 
EventScheduled(0) = False 
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End If 


" Check if the server is busy. 
If ServerBusy = 1 Then 


‘Server is busy, so put this customer at the end of the queue. 
NumiInQ = NumInQ + 1 


‘If the queue is now longer than it has been before, update MaxNumInQ and redimension arrays 
" appropriately. 
If NumInQ > MaxNumInQ Then 
MaxNumiInQ = NumInQ 


‘The "+1" in the next line is because QTimeArray is 0-based, so its elements are now 0 to 
MaxNumiInQ. 
ReDim Preserve QTimeArray(MaxNumInQ + 1) 


‘ TimeOfArrival is 1-based, with elements 1 to MaxNumInQ. 
ReDim Preserve TimeOfArrival(1 To MaxNumInQ) 
End If 


' Keep track of this customer's arrival time (for later stats). 
TimeOfArrival(NumInQ) = ClockTime 


Else 
"The customer can go directly into service, so update the status of the server to busy. 
ServerBusy = 1 


' Schedule a departure event for this server. 
EventScheduled(1) = True 
TimeOfNextEvent(1) = ClockTime - MeanServeTimeDivByK * Log(ProdKRand) 


End If 
End Sub 


Sub DepartureQ) 
Dim TimelInQ As Single, i As Integer 


"Update number of customers who have finished. 
NumServed = NumServed + 1 


' Check if any customers are waiting in queue. 
If NumInQ = 0 Then 


‘No one is in the queue, so make the server who just finished idle. 
ServerBusy = 0 
EventScheduled(1) = False 
Else 


" At least one person is in the queue, so take customer from front of queue into service. 
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NumInQ = NumInQ - 1 


"TimeInQ is the time this customer has been waiting in line. 
TimeInQ = ClockTime - TimeOfArrival(1) 


' Check if this is a new maximum time in queue. 
If TimeInQ > MaxTimeInQ Then 
MaxTimelInQ = TimeInQ 
End If 


"Update the total of all customer queue times so far. 
SumOfQTimes = SumOfQTimes + TimeInQ 


Schedule departure for this customer. 
TimeOfNextEvent(1) = ClockTime - MeanServeTimeDivByK * Log(ProdKRand) 


‘Move everyone else in line up one space. 
For i= 1 To NumInQ 
TimeOfArrival(i) = TimeOfArrival(i + 1) 
Next 
End If 
End Sub 


‘Sub Report0) 
Dim i As Integer, AvgTimeInQ As Single, AvgNumInQ As Single, AvgServerBusy As Single 


' Calculate averages. 
AvgTimeInQ = SumOfQTimes / NumServed 
AvgNumInQ = TotalTimeInQ / ClockTime 
AvgServerBusy = TotalTimeBusy / ClockTime 


* QTimeArray records, for each value from 0 to MaxNumInQ, the percentage of time that many 
customers were 
‘ waiting in the queue. 
For i = 0 To MaxNumInQ 
QTimeArray(i) = QTimeArray(i) / ClockTime 
Next 


"Enter simulate results in named ranges. 
Range("FinalTime") = ClockTime 
Range("NumServed") = NumServed 
Range("AvgTimeInQ") = AvgTimeInQ 

* Range("MaxTimeInQ") = MaxTimeInQ 
Range("AvgNumInQ") = AvgNumInQ 
Range("MaxNumInQ") = MaxNumInQ 
Range("AvgServerUtil") = AvgServerBusy 


‘Enter the queue length distribution from row 27 down, and name the two columns. 
With Range("A27") 
For i = 0 To MaxNumInQ 
.Offset(i, 0) =i 


.Offset(i, 1) = QTimeArray(i) 
Next 
Range(.Offset(0, 0), .Offset(MaxNumInQ, 0)).Name = "NumInQ" 
Range(.Offset(0, 1), .Offset(MaxNumInQ, 1)).Name = "PctOfTime" 
End With 


Range("A2").Select 
End Sub 


Sub UpdateChart() 
With Charts("QChart") 
.Visible = True 
Activate 
End With 
With ActiveChart 
With .SeriesCollection(1) 
.Values = Range("PctOfTime") 
-XValues = Range("NumInQ") 
End With 
Deselect 
End With 
End Sub 


Sub ViewChangelInputs() 
Call ClearOldResults © 
Call ViewReport 
End Sub 


Sub ViewReport(Q) 
With Worksheets("Report") 
. Visible = True 
Activate 
End With 
_Range("A2").Select 
End Sub 
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Appendix D 


Simulation Results 
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Stow Copier Averages 


Wy 
Runt 
Run 2 
Run 3 
Run 4 
Run 5 
Run 6 
Run 7 
Run 8 
Run 9 
Run 10 
Run 11 
Run 12 
Run 13 
Run 14 
Run 15 
Run 16 
Run 17 
Run 18 
Run 19 
Run 20 
Run 21 
Run 22 
Run 23 
Run 24 
Run 25 
Run 26 
Run 27 
Run 28 
Run 29 
Run 30 
Run 31 
Run 32 
Run 33 
Run 34 
Run 35 
Run 36 
Run 37 
Run 38 
Run 39 
Run 40 
Run 41 
Run 42 © 
Run 43 
Run 44 
Run 45 
Run 46 
Run 47 
Run 48 
Run 49 
Run 50 


0.10927 
0.16456 
0.03091 
0.02018 
0.05388 
0.04864 
0.06316 
0.16538 
0.06414 
0.30133 
0.04675 
0.01716 
0.01632 
0.04170 
0.01011 
0.01913 
0.04454 
0.20417 
0.07917 
0.91100 
1.07791 
0.42537 
0.03478 
0.07905 
0.17435 
0.08241 
0.13850 
0.11912 
0.16659 
0.43531 
0.14946 
0.44913 
0.04688 
0.12719 
0.39476 
0.07604 
0.18358 
0.37731 
0.67135 


0.34790 


0.09144 
0.55159 
0.02253 
0.03967 
0.04287 
0.15722 
0.65186 
0.01664 
0.10607 
0.40448 


0.50851 
0.60034 
0.14313 
0.07852 
0.30517 
0.25867 
0.18801 
0.56967 
0.30595 
1.22444 
0.24258 


* 0,06160 


0.06194 
0.14388 
0.04398 
0.07145 
0.23207 
0.72253 
0.29356 
3.51945 
3.48883 
1.95144 
0.29693 
0.31373 
0.69373 
0.37815 
0.75591 
0.42607 
0.65557 
1.49719 
0.53375 
1.89058 
0.16771 
0.49354 
1.92056 
0.45679 
0.68806 
1.55858 
3.16935 
1.05819 
0.41338 
2.02181 
0.21446 
0.16618 
0.24161 
0.51674 
2.21508 
0.07117 
0.35714 
1.44451 


0.11434 
0.05104 
0.35921 
0.01063 
0.25623 
0.07176 
0.04190 
0.09911 
0.05470 
0.20236 
0.12273 
0.01644 
0.24613 
0.08115 
0.29422 
0.04323 
0.08142 
0.05520 
0.26868 
0.11231 
0.57751 
0.06814 
0.18116 
0.18295 
0.45038 
0.24432 
0.03234 
0.76411 
0.12394 
1.03363 
0.04365 
0.07425 
0.21581 
0.22772 
0.14121 
0.17969 
0.78580 
0.19505 
0.01026 
0.09084 
0.12865 
0.28878 
0.11011 
0.12379 
0.13881 
0.09893 
0.50984 
0.09218 
0.01715 
0.28809 


0.49613 Average W,= 
0.23425 Average L, = 
1.12355 
0.04862 
1.15772 
0.23782 
0.21869 
0.34729 
0.21448 
0.65605 
0.32324 
0.05899 
1.10992 
0.32293 
1.04758 
0.14670 
0.41760 
0.21925 
1.15016 
0.47921 
2.78323 
0.21888 
0.54248 
0.76858 
1.48173 
1.06087 
0.12867 
3.19705 
0.75589 
3.95567 
0.27297 
0.31905 
1.22555 
1.16232 
0.44395 
0.68933 
2.71905 
0.74742 


~ 0.03390 


0.33656 
0.44078 
1.20993 
0.41508 
0.45177 
0.63914 
0.37081 
2.04580 
0.43785 
0.07032 
1.53465 


0.20055 
0.80102 


Fast Copier Averages 


W, 

Runi 
Run 2 
Run 3 
Run 4 
Run § 
Run 6 
Run7 
Run 8 
Run 9 
Run 10 
Run 11 
Run 12 
Run 13 
Run 14 
Run 15 
Run 16 
Run 17 
Run 18 
Run 19 
Run 20 
Run 21 
Run 22 
Run 23 
Run 24 
Run 25 
Run 26 
Run 27 
Run 28 
Run 29 
Run 30 
Run 31 
Run 32 
Run 33 
Run 34 
Run 35 
Run 36 
Run 37 
Run 38 
Run 39 
Run 40 
Run 41 
Run 42 
Run 43 
Run 44 
Run 45 
Run 46 
Run 47 
Run 48 
Run 49 
Run 50 


0.16522 
0.13201 
0.04372 
0.02195 
0.01193 
0.05944 
0.14870 
0.00621 
0.03050 
0.07030 
0.01497 
0.07449 
0.01731 
0.00629 
0.02977 
0.01510 
0.02951 
0.04862 
0.00509 
0.13592 
0.02036 
0.14392 
0.04557 
0.01441 
0.02689 
0.02470 
0.01629 
0.03730 
0.02351 
0.01063 
0.07228 
0.02899 
0.07347 
0.04713 
0.01518 
0.01361 
0.00924 
0.01632 
0.01872 
0.01248 
0.02905 
0.01285 
0.01141 
0.05441 
0.00999 
0.09387 
0.01463 
0.05637 
0.02227 
0.03851 


0.67717 
0.52072 
0.20845 
0.08367 
0.03721 
0.23190 
0.48656 
0.02343 
0.09118 
0.28705 
0.06099 
0.27182 
0.07681 
0.02254 
0.07318 
0.06749 
0.11189 
0.19437 
0.02329 
0.67674 
0.07768 
0.63696 
0.19706 
0.05882 
0.11062 
0.10410 
0.08179 
0.15623 
0.09129 
0.03970 
0.29630 
0.13787 
0.45476 
0.13813 
0.04490 
0.05502 
0.03200 
0.04468 
0.05931 
0.06478 
0.11359 
0.04996 
0.05132 
0.22992 
0.03312 
0.27498 
0.04960 
0.24002 
0.11487 
0.18161 


0.10561 
0.00541 
0.00838 
0.05808 
0.01744 
0.11641 
0.02950 
0.03575 
0.07305 
0.01799 
0.01810 
0.03163 
0.02978 
0.02474 
0.03375 
0.07908 
0.00127 
0.03146 
0.08552 
0.01336 
0.00732 
0.03880 
0.04772 
0.01123 
0.00153 
0.01957 
0,02922 
0.02355 
0.03911 
0.15835 
0.02046 
0.04936 
0.01360 
0.03550 
0.03414 
0,02918 
0.01492 
0,09993 
0.02777 
0.07698 
0.00724 
0.03041 
0.13834 
0.00983 
0.11393 
0.01693 
0.07004 
0.05218 
0.01722 
0.02611 


0.44728 Average W,. 
0.01846 Average L, = 
0.02305 
0.20816 
0.08226 
0.92124 
0.11186 
0.13533 
0.31083 
0.05885 
0.06762 
0.10642 
0.02346 
0.03843 
0.19256 
0.28537 
0.00649 
0.10932 
0.16593 
0.06215 
0.02680 
0.09256 
0.22064 
0.05152 
0.00650 
0.06733 
0.06145 
0.01329 
0.14490 
0.80359 
0.11354 
0.04148 
0.06785 
0.06041 
0.17413 
0.10268 
0.01783 
0.30203 
0.09101 
0.54003 
0.03081 
0.13233 
0.89531 
0.03708 
0.53144 
0.06160 
0.26666 
0.24669 
0.08720 
0.13270. 


0.04158 
0.16644 
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